This notebook tries to illustrate the concept of a rare event, the logic behind employing an umbrella potential, and Metadynamics (both the plain and the well-tempered version). Everything is implemented trying to keep things as simple as possible. The notebook relies on numpy, matplotlib, scipy, and tqdm.
import numpy as np
%matplotlib inline
import matplotlib as mpl
import matplotlib.pyplot as plt
from tqdm import tqdm
# Set a random seed for reproducibility
np.random.seed(2021)
# Resolution for the various images
mpl.rc('figure', dpi=200)
Everything will be performed on a small toy system. The system we consider is 1D. We use the following functional form for the potential energy:
\begin{equation} U(x) = ax^4 - 6x^2 + x \end{equation}with $a=0.5$, $b=-6$, and $c=1$. It follows that we have a very simple expression for the force:
\begin{equation} f(x) = -\frac{dU(x)}{dx} = -2x^3 + 12x - 1 \end{equation}def potential_energy(x):
'''Compute the value of the potential energy'''
return 0.5*x**4 - 6*x**2 + x
def force(x):
'''Compute the force from the potential energy'''
return -2*x**3 + 12*x - 1
The potential energy has the following look:
# Evaluate the potential energy over a grid spanning from -4 to +4 (100 bins)
xx = np.linspace(-4,4,100)
U = potential_energy(xx)
def plot_U(xx, U):
'''Plot the potential energy'''
plt.plot(xx, U, lw=3.0, c='k')
plt.xlim(-5,5)
plt.ylim(-25,25)
plt.tight_layout()
plot_U(xx, U)
Then we define a simple System class. It's nothing more than a container with the initial coordinates and an attribute traj that will store the trajectory
class System():
def __init__(self, x0, p0, m, kbT):
self.x0 = x0
self.p0 = p0
self.m = m
self.kbT = kbT
# Also initialize a list that will
# store the trajectory.
self.traj = [x0]
We can now initialize our system. We initialize it in the basin on the right (x0 = 2.4), with no initial velocity (p0 = 0.0), with mass equal to one (m = 1), and with a small $k_bT$ (kbT = 1e-2)
# Settings for the `System`
x0 = 2.4
p0 = 0.0
m = 1.0
kbT = 1e-2
# Initialize the system
system = System(x0=x0, p0=p0, m=m, kbT=kbT)
We can take a look at it
def plot_point_on_U(x, xx, U, shift=2.3):
'''Plot the system over U'''
plot_U(xx,U)
plt.plot(x,potential_energy(x0)+shift, 'o', markersize=20, markeredgecolor='k', color='darkorange')
plot_point_on_U(system.x0, xx, U)
Here we define a simple class that represents a molecular dynamics engine. We simulate our system in the NVT ensemble, by means of the Langevin equations of motions (Langevin dynamics). We integrate the equations of motions using the BAOAB method.
\begin{align} p(t+\frac{1}{2}\delta t) &= p(t) + \frac{1}{2}\delta t f(f) &B\\ x(t + \frac{1}{2}\delta t) &= x(t) + \frac{1}{2}\delta t p(t + \frac{1}{2}\delta t)/m &A \\ p'(t+\frac{1}{2}\delta t) &= \exp(-\gamma \delta t)p(t + \frac{1}{2}\delta t) + \sqrt{1-\exp(-2\gamma\delta t)}\sqrt{mk_bT}G &O \\ x(t+\delta t) &= x(t + \frac{1}{2}\delta t) + \frac{1}{2}\delta t p'(t + \frac{1}{2}\delta t)/m &A \\ p(t + \delta t) &= p'(t+\frac{1}{2}\delta t) + \frac{1}{2}\delta t f(f + \delta t) &B \end{align}Here $G$ is a Gaussian random variable with zero mean and unit variance, and $\gamma$ is the frictional coefficient.
The class is initialized with some MD settings like the gamma constant of the Langevin dynamics, the time step, the number of steps to do, and the function that evaluates the force. It has a method run that runs the simulation and stores the coordinates in the given System.
class MDEngine():
def __init__(self, langevin_gamma, dt, n_steps, force):
self.langevin_gamma = langevin_gamma
self.dt = dt
self.n_steps = n_steps
self.force = force
def run(self, system):
'''Run the simulation for the provided system'''
from tqdm import tqdm
# Initialize the variables and compute
# the initial force.
self._x = system.x0
self._p = system.p0
self._f = self.force(system.x0)
self._m = system.m
self._kbT = system.kbT
# Run the dynamics
for step in tqdm(range(self.n_steps)):
self.langevin_step()
# Also update the system trajectory
system.traj.append(self._x)
def langevin_step(self):
'''Langevin step using BAOAB integration (see Allen and Tildesley)'''
self._p = self._p + 0.5*self.dt*self._f # B
self._x = self._x + 0.5*self.dt*self._p/self._m # A
self._p_r = np.exp(-self.langevin_gamma)*self._p + np.sqrt(1-np.exp(-2*self.langevin_gamma*self.dt))*np.sqrt(self._m*self._kbT)*np.random.normal() # O
self._x = self._x + 0.5*self._p_r/self._m # A
self._f = self.force(self._x)
self._p = self._p_r + 0.5*self.dt*self._f # B
Now we can initialize the system and the engine. We first consider the case of a high barrier. "High" here refers to the size of the barrier compared with $k_bT$, i.e., we set $k_bT$ to a small value.
# Settings for the `System`
x0 = 2.4
p0 = 0.0
m = 1.0
kbT = 0.5
# MD Engine settings
gamma = 1.0
dt = 5e-3
n_steps = 100000
force = force
# Initialize the system and the engine
system = System(x0=x0, p0=p0, m=m, kbT=kbT)
md_engine = MDEngine(langevin_gamma=gamma, dt=dt, n_steps=n_steps, force=force)
# Run the simulation
md_engine.run(system)
100%|██████████| 100000/100000 [00:01<00:00, 76304.82it/s]
If we plot the trajectory, we see that the system visits only the basin on the right
def plot_traj(traj):
plt.plot(traj, c='darkorange')
plt.xlabel('Time')
plt.ylabel('s')
plot_traj(system.traj)
We can also visualize the trajectory directly. This may take some time, and more so if the trajectory gets longer.
def visualize_dyn(traj, stride=100, shift=2.3):
from matplotlib import animation
from IPython.display import HTML
fig, ax = plt.subplots(1, 1)
ax.set_xlim(-4,4)
ax.set_ylim(-25,25)
ax.plot(xx, U)
def animate(frame):
ax.cla()
ax.plot(xx, U, lw=3, c='k')
ax.plot(traj[::stride][frame], potential_energy(traj[::stride][frame])+shift, 'o',
c='darkorange', markersize=20, markeredgecolor='k')
ani = animation.FuncAnimation(fig, animate, frames=len(traj[::stride]))
return ani
ani = visualize_dyn(system.traj, stride=500)
from IPython.display import HTML
HTML(ani.to_jshtml())